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Understanding the interaction mechanisms between large-scale magnetohydrodynamic instabili¬ 
ties and small-scale drift-wave microturbulence is essential for predicting and optimizing the per¬ 
formance of magnetic confinement based fusion energy experiments. We report progress on under¬ 
standing these interactions using both analytic theory and numerical simulations performed with 
the BOUT-h+ [B. Dudson et ah, Comput. Phys. Comm. 180 , 1467 (2009)] framework. This 
work focuses upon the dynamics of the ion temperature gradient instability in the presence of a 
background static magnetic island, using a weakly electromagnetic two-dimensional hve-held fluid 
model. It is found that the island width must exceed a threshold size (comparable to the turbulent 
correlation length in the no-island limit) to signihcantly impact the turbulence dynamics, with the 
primary impact being an increase in turbulent fluctuation and heat flux amplitudes. The turbulent 
radial ion energy flux is shown to localize near the X-point, but does so asymmetrically in the 
poloidal dimension. An effective turbulent resistivity which acts upon the island outer layer is also 
calculated, and shown to always be significantly (lOx - lOOx) greater than the collisional resistivity 
used in the simulations. 


I. INTRODUCTION 

One of the most crucial factors in determining the size 
of magnetic-confinement based fusion energy (MFE) re¬ 
actors is the level of plasma confinement achieved [1]. In 
MFE-relevant plasmas, the dominant physical processes 
which determine this confinement are a variety of insta¬ 
bilities driven by the inherent temperature, density, and 
current gradients of the confined plasma. Two of the 
most important such instabilities in tokamaks are the ion 
temperature gradient (ITG) mode [2] and tearing modes 
[3]. As its name suggests, the ITG instability is driven by 
the radial gradient (i.e. across nested magnetic flux sur¬ 
faces) of the equilibrium ion temperature profile, which 
must exceed a critical value (the so-called critical gradi¬ 
ent) for onset of the instability. Once the critical gra¬ 
dient is exceeded, the ITG instability leads to a broad 
spectrum of ion gyroradius {pi) scale fluctuations, whose 
collective nonlinear behavior drives rapid cross-field par¬ 
ticle, energy, and momentum fluxes that greatly exceed 
collisional transport levels [4], and often determine the 
overall confinement level achieved. Since these fluctua¬ 
tions generally have correlation lengths on the order of 
1 — lOpi, which is much smaller than the minor radius a 
in MEE-relevant devices, the saturated fluctuations and 
their dynamics are often referred to as microturbulence. 

Magnetic islands can appear through magnetohydro¬ 
dynamic (MHD) instabilities [5, 6], such as the tear¬ 
ing mode (TM) driven by equilibrium gradients, or be 
externally imposed with resonant magnetic perturba¬ 
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tions (RMP). These islands can reduce the confinement 
achieved, similar to the ITG instability, or even termi¬ 
nate confinement by leading to a disruption. It is theo¬ 
retically predicted and experimentally observed that the 
presence of a magnetic island flattens the electron and ion 
temperature and density. This flattening occurs both in¬ 
side the island through rapid parallel equilibration, both 
through parallel sound waves [7] and rapid parallel heat 
conduction [8], and in time averaged profiles by enhanc¬ 
ing the effective cross-field transport. The classic picture 
is that the radial heat flux increases at the X-point of 
the magnetic island due to the fast parallel dynamics in 
the highly anisotropic 3D island structure [8]. Indeed, 
because of the direct connection of the magnetic flux 
surface at the separatrix of the island, all density and 
temperature structures on one side of the island can be 
rapidly transported to the other side of the island on 
shorter timescales than the timescale of the background 
cross-field transport. 

Given the importance of both classes of instabilities in 
determining the level of confinement achieved in magnet¬ 
ically confined plasmas, it is important to develop mod¬ 
els of their dynamics which can be used to accurately 
interpret existing observations and confidently design fu¬ 
ture experiments and reactors. While great progress has 
been made in recent years in understanding the nonlin¬ 
ear dynamics and scalings of each class of instabilities in 
isolation (often through use of massively parallel compu¬ 
tation), these instabilities frequently coexist in current 
high-performance tokamak discharges, and are expected 
to continue to do so in many future reactor scenarios. 
We must therefore also understand how these instabili¬ 
ties couple linearly and nonlinearly if we are to be able 
to predict the overall confinement and performance of 
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such discharges. For example, it is now well-known that 
microturbulence (and ITG microturbulence in particu¬ 
lar) generally saturates via the nonlinear generation of 
axisymmetric zonal flows [9] which in turn shear apart 
the ITG eddies, forming a self-regulating system of tur¬ 
bulence and zonal flows. How this process is altered by 
the presence of magnetic islands (which generally have 
toroidal mode numbers n = 1 — 2 and poloidal mode 
numbers m = 2 — 5) and their associated flow fields re¬ 
mains an open question. Alternatively, the saturation of 
the neoclassical tearing mode depends sensitively on the 
level of profile flattening that occurs in and near the is¬ 
land [10, 11]. As this flattening depends sensitively on 
the ratio of parallel and perpendicular (microturbulence- 
dominated) transport, it is clear that the saturated island 
width will depend self-consistently on the level of micro¬ 
turbulence present (which also depends on the amount 
of profile flattening that occurs). 

These two instabilities (ITG and TM) are well studied 
in the literature, with many different research teams find¬ 
ing a wide array of complex nonlinear couplings between 
microturbulence and tearing modes. For instance, McDe- 
vitt and Diamond [12] demonstrated that drift-wave mi¬ 
croturbulence could nonlinearly excite an island via the 
Reynolds stress in a fashion analogous to the formation of 
zonal flows, while Sen et al [13] used a similar approach to 
demonstrate that electron temperature gradient (ETG) 
modes could nonlinearly damp the island via an effective 
turbulent resistivity. More recently, a wide range of re¬ 
searchers have utilized both fluid [14-26] and gyrokinetic 
[27-37] approaches to investigate couplings of various mi¬ 
croturbulence instabilities to either statically imposed is¬ 
lands or dynamically evolving tearing modes. Although 
there is no simple and clear consensus yet, it is clear that 
the presence of the island can lead to a poloidal localiza¬ 
tion of linear microturbulence eigenmodes, and that the 
microturbulence can both damp and destabilize islands 
depending on the type of microturbulence instability, and 
the width and rotation speed of the island. The impacts 
of the turbulence on the island evolution are generally 
represented by inclusion of turbulent polarization current 
terms in generalized Rutherford equations [18, 38]. 

Building upon these results, we present in this paper 
the findings of a study examining the responses of weakly 
electromagnetic ITG turbulence to statically imposed is¬ 
lands of varying width, performed using a simple fluid 
model in slab geometry. We find that above a critical 
island width, the magnitude of both the microturbulence 
fluctuations and associated radial ion heat flux Qi in¬ 
crease approximately linearly with island width, both 
near and far from the island resonant surface, for val¬ 
ues of the the ion temperature gradient both close to, 
and well above the linear critical gradient. The critical 
island width is found to be approximately equal to the 
radial correlation length of the turbulence in the absence 
of an island. Perhaps not surprisingly, for islands smaller 
than the critical size, the microturbulence is found to be 
effectively insensitive to the presence of the island. We 


also observe the turbulent Qi to localize near the island 
X-point as expected, but somewhat surprisingly, this lo¬ 
calization is not seen to be symmetric in the poloidal 
direction. To begin assessing how the turbulence may be 
expected to back-react upon the island, we present calcu¬ 
lations of an effective turbulent resistivity which operates 
on the “outer solution” of the island structure. We show 
that for the parameters considered, while this resistivity 
is always much larger than the collisional value, it max¬ 
imizes for small island size and decreases rapidly with 
increasing island width. 

The paper is organized as follows. Section II details the 
fluid model used to study the effects of a static magnetic 
island on the ITG microturbulence, and includes details 
of the numerical implementation in BOUT+-h (IIB) and 
verification of the implementation using analytic growth 
rate calculations (IIG) and a global energy balance anal¬ 
ysis suitable for nonlinear simulations (HD). Sec. Ill dis¬ 
cusses the results of numerical simulations quantifying 
the impact of statically imposed islands of varying width 
on ITG turbulence, and Sec. IV discusses initial work ex¬ 
ploring the expected feedback of the turbulence on the 
island. Gonclusions and future research directions are 
discussed in Sec.V 


II. DETAILS OF THE FIVE-FIELD FLUID 
MODEL 

This section describes the 2D electromagnetic five-field 
model used in this study, and its implementation in the 
BOUT-j—|- framework [39]. Two verification methods are 
presented to test the implementation. The first is a com¬ 
parison of calculated linear growth rates against analytic 
derivations, detailed in Sec. IIG. The second is a global 
energy balance analysis approach, detailed in Sec. IID. 

A. Model overview 

In order to investigate the self-consistent couplings of 
ITG turbulence and tearing modes in a simple, tractable 
form, we employ a simple five-field two-fluid model [21, 
40, 41] in a Gartesian slab geometry [21, 41-44] simi¬ 
lar to that used in previous studies, considered together 
only in Ishizawa et al. [21, 41]. The simulation radial 
and poloidal domains span the range —a<x<a and 
—6/2 < y < 6/2 respectively, with Dirichlet boundary 
conditions on all fluctuating fields at x = ±a and periodic 
boundary conditions at ^ = ±6/2. Using the standard 
definitions oi E = —V(j) and B = Bq{x)z ± Wip x z, the 
model describes the evolution of fluctuations in the elec¬ 
tron density n = Snjn^ vorticity D = Vic/., ion parallel 
velocity Vi^z = ion temperature Ti = 6Ti/T^ and 

magnetic flux ^ = {1 jp){cs jc)e5^pjT^ where 0 = e5(p/T 
is the normalized electrostatic potential fluctuation, and 
the normalized axial current is given hy jz = 

Here, n and T are constant normalizing factors taken 
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to be equal to the values of the equilibrium density and 
ion temperature profiles at x = 0; = ^JT/Mi and 

ps = Cs/O^cii with Oci = eBo(0)/M^c. In this model, 
we can chose a form of the equilibrium magnetic flux 
that allows us to specify both an equilibrium poloidal 
magnetic field By^Q{x) = —dx'ipo ^.nd a static magnetic 
island Bx,isiand{y) = dy^pisiand- These relations come 
from B = V'0(x,^) x z and the current is given by 
jzo = z'VxB = — To further simplify the problem, 
in the rest of this paper we restrict our attention to the 
case for which the equilibrium density no is constant in 
radius and the equilibrium electrostatic potential and 
parallel ion flow Vi^zO are equal to zero. The equilibrium 
ion temperature profile is taken to have a simple linear 
dependence Tio{x) = T{\—xjLTi)^ where Lti is constant 
across the radial domain. The equilibrium magnetic flux 
is taken to have the form y) = '0o(^)+'0isZand(^) = 

-x^l{2Ls) + V’iCos(27r2//6), with V’i = 
and the associated current is jzo{x^y) = —where 
= ^xx ^^yy' constituent equations of the model 
are the electron density, the perpendicular momentum 
associated with the quasi-neutrality, the parallel momen¬ 
tum, the ion temperature conservation, and Ohm’s law 
as follows: 

dtfi = 


dt^ = 


+ 

dfXi^z — 
dtfi = 


Pdd = 


V^dy 


- h) + V|| {jz - Vi^z) + Viii^o 


(^)t 


D I h 


,y ■ 

Viiiz + V||j^o 

V*Bdy{{l + T)h + TT,} + TV^.dyh 


rV 


_L • 


Diy\Q, 


-V|| {p.+T)h + Tf^ 

Vs {Vi,z)y + D±V]_Vi^z 


Xv%{T - l)dy 


T + n + Ti 


Vs (Ti) + Z)_lVj_T, 


-Vii 


- n -r]j^ 


( 1 ) 

( 2 ) 

( 3 ) 

( 4 ) 

( 5 ) 

(6) 


where = a/L b and = a/LTi are the driven im¬ 
posed gradients of the magnetic field and the ion tem¬ 
perature. In these equations x and y are normalized to 
ps and t to a/cs- Using the standard Poisson bracket 
notation {f,g} = djjdxdgjdy — djjdydgjdx, we write 

the total time derivative as djjdt = dfjdt + |^,/| to 
include the E x B convection term, and the total paral¬ 
lel derivative V||/ = -(3 + Asiand + i>, /| includes 

derivatives along the total magnetic field (the equilib¬ 
rium field, any imposed island, and the magnetic field 
fluctuations). We also note that the use of the notation 

/| is associated only with the magnetic 
fluctuations part. The coefficient T = 5/3 is the ratio 


of specific heats and r is the ratio between the ion and 
electron reference temperatures. The coefficient A, equal 
to 0 or 1, is used to investigate the impact of different 
closure models in the Ti equation. We have included A 
in order to investigate the effects of the absence in some 
models found in the literature of the divergence of the 
ion E X B and diamagnetic velocity terms. This term, 
describing finite Larmor-radius (FLR) effects, is retained 
(or not) in the fluid closure model used for the ion tem¬ 
perature equation. Moreover, this model contains other 


FLR effects (i.e., r'V±_ • 


V_L0, ) 


[45-47] which 


are not always considered in the literature but included 
here in order to take into account the ion temperature ef¬ 
fects which come from the convection of the ion pressure 
by the E x B drift. This term appears in the vortic- 
ity equation due to the divergence of the ion polariza¬ 
tion current. When the inhomogeneity of the magnetic 
field is kept, the models including this FLR term such 
as in Refs. [40, 45, 46, 48-51] are not Hamiltonian [52] 
in the ideal limit {D = y = Ug = 0). However, the 
ideal version of our model is very close to a Hamilto¬ 
nian structure because the terms required to recover the 
Hamiltonian property involve, for example, higher order 
spatial derivatives of the magnetic field which are can¬ 
celed by the assumption of a linear background magnetic 
gradient, or higher order spatial derivatives of the vor- 
ticity which are assume negligible. Moreover, the energy 
conservation of our simulations is systematically verified 
as well as the error between the numerical and theoret¬ 
ical computation of the energy balance. We found that 
the inclusion of the FLR term (with A = 1) changes the 
growth rate less than 10% in the linear regime, but this 
FLR term in the ion temperature equation and the one 
in the vorticity equation are kept in the following be¬ 
cause they contribute to the numerical stability of our 
simulations for high ky modes. 

This set of equations becomes an electrostatic two- 
dimensional four-field fluid model with gradient driven 
turbulence when /3 = 0. Moreover, when (3 is relatively 
small (i.e., ^ 10“^) the differences between the five-field 
and the four-field electrostatic model are not significant. 
The five-field model is preferred for consistency with fu¬ 
ture work on dynamic islands. The additional dissipative 
viscosity Ug is included to prevent quasilinear relaxation 
and maintain input background profiles such as the back¬ 
ground ion temperature gradient. Without this forcing, 
the ion temperature gradient cannot be maintained in the 
presence of a magnetic island and the ITG which drives 
the turbulence becomes stable. This choice is needed to 
maintain on average the same level of turbulence due to 
the ITG. We show below that this assumption is still 
compatible with the flattening effect of the island on the 
temperature and density profiles. The numerical dissipa¬ 
tion terms proportional to utilize small values of 

to damp grid-scale fluctuations without significantly 
impacting the large-scale dynamics of interest. Note that 
relative to previous tearing mode studies which have uti¬ 
lized anomalous transport coefficients, these simulations 
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self-consistently determine those coefficients via the tur¬ 
bulent fluxes e.g. Qi = (fiVExs). 


B. Numerical implementation 

The BOUT- 1 -+ [39] framework is used for the numer¬ 
ical computation of the five-field fluid model given by 
Eqs. (l)-( 6 ). Even though BOUT++ can deal with a 
helical axisymmetric magnetic field, as a first step our 
work focuses on a slab geometry choosing the box ratio 
of the simulations in order to fix a specific regime of the 
instabilities. The use of BOUT-I-+ is an asset since the 
dynamical equations are written in a small script with 
object oriented functions. Also, all multi-dimensional 
scans and post-treatment analyses have been simplified 
by the use of the OMEIT [53] framework and has mo¬ 
tivated the creation and development of the BOUT- 1 — 1 - 
interface module in the OMEIT framework as explained 
in Ref. [54]. As a summary, OMEIT is a software (with 
Graphical User Interface) for integrated studies which 
can quickly connect many codes and experimental data 
in a workflow and which can perform data analysis. Eor 
instance, OMEIT was used to handle remote execution, 
data transfer, and analysis of the simulations presented 
here, which were performed mainly on the TSCC com¬ 
puter at the San Diego Supercomputer Center. 


C. Linear dispersion relation 

In order to verify the accuracy of our numerical sim¬ 
ulations against the target fluid model, the first step is 
to compare the computed linear growth rates against an¬ 
alytic calculations in the limit of no background static 
island and no static magnetic shear (i.e., 'ipeq = 0 ). 
The linearization of the fluid equations by f{x^y^z^t) = 
f exjp{ikxX + ikyU — iut) yields 


— ujfl 

1 

II 

(7) 

ujk\^ 

= -v*Bky{{l + T)h +rfi) 



—TVB.kyk\(j), 

(8) 


= 0, 

(9) 

-Ujfi 

II 

1 



-\-v^kyX(T — 1) ^0/t + h -h Ti^ , 

(10) 


= iyk^ip. 

(11) 


After some analytic computation the dispersion relation 
is reduced to 


with 

As = fikl, (13) 

A 4 = k'i{i7]k'i+ISC bti), (14) 

As = -l3flBii'i- + T + r2)flB-TflTi) + k‘]_iriCBTi 

-l3klnB{r2nB - r(r 2 - mn), (is) 

A2 = -ir]kj_{{l + T+ T2)Q.% - TQ.TiO.B - kj_CBTi) 

-l3Ti}Ti^%{T2k‘i + 1 ), (16) 

Ai = kj_iri(r2k'j_ + l)TQ,Ti^%, (17) 

Ao = 0 , (18) 

where fts = ^2 = A(r — 1 ) and 

CsTi = (r 2 — + rflTi- 

This linear dispersion relation is obviously driven by 
the imposed gradients of the ion temperature v^- and 
magnetic field v^. If the latter are dropped, the disper¬ 
sion relation becomes D{uj) = k^uj'^{Puj -h irjk‘]_) and the 
roots are found to be stable (imaginary part of cj, i.e. 
the growth rate, is 0 and the propagation frequency is 
equal to 0 or —k\r]l^'). A scan of the ion temperature 
and magnetic field gradients has shown that the linear 
growth rate is unstable beyond a threshold dependent on 
both gradients. We typically perform our simulations in 
the unstable regime with = 0.1 and =0.3, highly 
relevant to experimental values. 

A comparison of the computed linear growth rates 
against the analytic linear dispersion relation (Eq. (12)) 
is shown in Eig. 1, demonstrating excellent agreement in 
both growth rate 7 and real frequency cj for this base 
case; similar agreement is obtained at other parameters. 
Error estimates on the simulation results are calculated 
via the standard deviation of the instantaneous frequency 
and growth rate of the time-averaging window used (gen¬ 
erally on the order of lOO’s of a/cg). Eor wave num¬ 
bers with sufficiently small or zero growth rate, a well- 
converged value is often not found, leading to the large 
error bars shown at the highest values of ky plotted and 
at ky = 0 . 




(a) Growth rate 7 (b) Propagation uj 

FIG. 1: Verification of the linear dispersion relation compu¬ 
tation between the theory (blue curves) given by the Eq. (12) 
and the BOUT-h-h simulation (dashed-black curves) of the 
Eqs. (l)-(6). The error bars (x-markers) of the simulation 
are plotted. The error bar becomes artificially big when the 
growth rate is stable or very small in the unstable regime. 


D{uS^ — -j- A4.(jj^ -|- A^uj^ -]- A2(jo‘^ -|- A\uj -|- Ao,(12) 
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D. Energy evolution and energy balance 


The next step of the verification process is to quan¬ 
tify the global energy balance of the simulations, which 
provides an error measure when the analytic calcula¬ 
tion of the linear dispersion relation cannot be used, due 
to inclusion of a non-uniform background magnetic field 
and/or nonlinear effects. The total energy of the model 
given by Eqs. (l)-( 6 ) is 


"= 2 


^ J (fx (1 + r) (^n 2 + |V_L(^p + /3|V_Lt^p) 


+v^ H- 


(19) 


where the first term corresponds to the internal poten¬ 
tial energy of the ions and electrons, the second term to 
the perpendicular kinetic energy, the third term to the 
magnetic fluctuation energy, the fourth term to the par¬ 
allel kinetic energy and the final term to the ion thermal 
energy. After some analytic computations the time evo¬ 
lution of the total energy is 


dfS = Ep -h Em + Eq + Ej — — Dp — (20) 

by specifying all components of the energy balance where 
the individual terms are given by 


Ep 

Em 

Eq 

Ej 

Dr, 


Dd 


D. 


XtvbP, 


(21) 

t(1 + t)v%M, 


(22) 

{r(l + r) A} ^ 4 , 

Q, 

(23) 

(1 + r) ( (n - ^) V\\jzo{x, 2 /))^ ^ 

(24) 

(l + r)r? 

D^{{1 + t) (|Vin|2 + |Vi0|2^ 

) 

(25) 

+ \v jfA 

/ x,y 

(26) 

Vs [[(1 + r) {h)l^y 



+ ('^hz)x,y + r - 1 ) ’ 


(27) 


with the positive driven terms related to: a pressure 
gradient flux P = {Tidyfi) , the radial particle flux 

\ / x,y 

M = {nvx)j.y, the radial ion heat flux Q = (TiV^) 

due to the E x B drift [55] where Vx = —9^0, and the 
negative dissipative terms related to: the resistivity 77 , 
the numerical diffusion D± and the axisymmetric sink 
terms proportional to which are used for the conser¬ 
vation of the background radial profiles. The coupling of 
the fluctuations to background axial and island current 
gradients Ej is added in order to verify its contribution 


in the energy balance. However, our choice of magnetic 
shear Byo{x) does not contribute to Ej since its second 
derivative vanishes, and so for the cases considered here, 
Ej contains only the contribution of the imposed island 
current gradient. In Fig. 2 and 3 below it is shown to be 
negligible for all island widths relative to the dominant 
source and sink terms. 

Each of these terms can be calculated in the simula¬ 
tions and the energy conservation properties of a given 
simulation can thereby be quantified in detail. The terms 
Ef are in general positive drives {Ej can be positive or 
negative), and the terms Df are dissipative contributions. 
As mentioned previously, we note that for the equations 
for Ep and Eq^ the presence of the FLR term A = 1 
stabilizes the turbulence. In fact, the magnetic gradient 
part of the energy balance related to the heat flux Eq 
is divided by 2 (when r = 1 ) and Fig. 3 shows a nega¬ 
tive contribution of Ep due to the presence of the FLR 
term. By calculating each of these terms from the simu¬ 
lation output, a numerical error measure e can be readily 
defined as 


^tEfium dfS 


For the nonlinear simulations we are using 125 radial 
grid points and 128 poloidal grid points with poloidal pe¬ 
riodicity and radial Dirichlet boundary conditions. The 
size of the box is (2a, 5) = (65ps,65ps) so the discretiza¬ 
tion is (dx^dy) ^ (0.52^^, 0.50ps). There are no signifi¬ 
cant effects of increasing the number of grid points by a 
factor 2 in each direction. With these parameters, the mi¬ 
croturbulence is well described and the saturation state 
is reached after a few hundreds of time steps (a/cg) as 
shown below. 

Fig. 2 shows the time evolution of the accuracy of the 
energy balance between the analytic prediction given by 
Eq. (20) and the numerical computation. All components 
of the analytic prediction of the energy balance are shown 
in different colors. We can clearly see that the dominant 
component of the energy balance is the one related to the 
radial ion heat flux, as expected in our ITG driven model. 
The numerical computation of the energy balance given 
by the time derivative of the numerical total energy is 
almost exactly identical to the analytic prediction. The 
relative error between the analytic and numerical com¬ 
putation of the energy balance e given by Eq. (28) is 
negligible with respect to the dominant term Eq . We do 
not show here but, as expected, the part of the energy re¬ 
lated the ion temperature fluctuations is dominant with 
respect to the other terms in Eq. (19). This explains the 
dominance of Eq with respect to the other components 
Ep^ Em and Ej. The values of the dissipative terms 
{D = 0.02, Us = 0.5, 77 = 10“^, and (3 = 10“^) act on the 
saturation of the numerical simulations. Indeed in this 
nonlinear simulation, with the chosen parameters, the 
saturation is obtained after approximately 400 a/cs time 
units since the energy balance fluctuates around zero for 
a static island of width Wisiand = 15ps. 
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In Fig. 3, time-averages (over the steady-phases) of 
normalized energies and energy balances dfS with vary¬ 
ing island width and the ion temperature gradient 
are shown, demonstrating the maintenance of good con- 


0.06 



^(a/cj 


FIG. 2: Time evolution of the energy balance in arbitrary 
units with a static background magnetic island of width 
Wisiand = 15ps. Each term of Eq.(20) is represented in color 
as well as their sum dtS (plain black curve) and the numerical 
value dfEnum (dashed-black curve). The relative difference e 
defined by Eq. (28) between the two later is negligible (bold 
black curve). 



"^islandiPs ) 


EIG. 3: Scaling in arbitrary units of time-averaged terms in 
the energy balance analysis as a function of v^. = ajLri over 
the time range t G [500,3000] a/cs in the presence of islands 
of widths Wisiand = {0,1,2,5,10,15}, and magnetic island 
size at fixed v^. = {0.2,0.25,0.3}. The sum of the positive 
source (respectively negative dissipative) terms of the energy 
balance given by Eq. (20) are noted S (respectively D). The 
error bars correspond to the standard deviation. 


servation for the simulation results discussed in the fol¬ 
lowing sections. All curves correspond to an ion temper¬ 
ature gradient of v^, = 0.3 except the yellow and cyan 
curves which are the dominant term of the energy bal¬ 
ance Eq for v!^, = 0.2 and 0.25. The numerical error 
is shown by the difference dt£ = S — D between the 
source S = Ep + Em + Eq + Ej and dissipative terms 
D = Dp Djy which is near 0. In particular, the 

value of the standard deviation of the theoretical predic¬ 
tion of the energy balance (std(dtf)) averaged over the 
island width is ~ 5 x 10 “^ which is much smaller than 
the standard deviation of the dominant term (std(F^Q)) 
averaged over the island width ^ 10“^. Moreover, the 
theoretical prediction of the energy balance dfS is not 
0 but is around 10“^, which is much smaller than the 
dominant term Eq ^ 0.03. Saturated energies increase 
with island width as discussed in Sec. Ill, however the 
dependence of the normalized energy and energy balance 
on island width is weak, and the numerical error is low. 


III. QUANTIFYING THE IMPACT OF STATIC 
ISLAND ON ITG TURBULENCE 

In our simulations we force the system to maintain 
on average the same level of turbulence due to the ITG 
by keeping an average background ion temperature 
gradient. The electrons are assumed to be correlated 
with the ions and are not the focus of this work. The 
time evolution of the magnetic island is known to be 
at least one order of magnitude slower than the time 
evolution of the microturbulence driven by the ITG. 
An electrostatic model {/3 = 0) allowing the exact force 
balance of the fluctuations in the Ohm’s law between 
the parallel classical resistive current and the parallel 
gradient of the drifts (i.e., the electric and diamagnetic 
drifts) is not used in this paper. Instead, weakly 
electromagnetic (small ^ = 10 “^) cases are considered, 
close to the electrostatic cases but, preparing for the 
full electromagnetic self-consistency simulations and 
allowing comparisons with other published electrostatic 
results. The advantage of full electromagnetic sim¬ 
ulations is that the effects of the turbulence on the 
slow dynamical evolution of the magnetic island will 
be self-consistently included without using the artifi¬ 
cially large numerical perpendicular transport coefficient. 

In Fig. 4 the saturated energy F^sat averaged on the 
time range t G [500; 3000] a/cg is plotted against the 
size of the static magnetic island as well as its stan¬ 
dard deviation. We observe that beyond a threshold 
of the island size (around bps) the saturated energy in¬ 
creases significantly. Further analysis has shown that it 
is due to an increase of the ion temperature fluctuations 
and not to the island magnetic energy. Indeed, as men¬ 
tioned previously, the dominant part of the energy given 
by Eq. (19) is the fluctuations of the ion temperature 
1/2/ rT//(r — 1). The values of 7} fluctuate around 
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FIG. 4: Energy saturation in arbitrary units for different is¬ 
land sizes averaged over the time range t G [500; 3000] a/cs 
for an ion temperature gradient = 0.3 (black curve), 0.25 
(blue curve) and 0.2 (red curve). The error bars correspond 
to the standard deviation. 


0.4 with no island and around 0.7 with Wisiand = 

The amplitude of the fluctuations are approximately 2 
times bigger (i.e., ^ ±0.2) for Wisiand = '^^Ps with re¬ 
spect to the one with no island. Fig. 5 shows time slices 
at times t = {1506,1516,1526,1536,1546}a/cs of the ion 
temperature T^o ± Ti without (Wisiand = Ops) and with 
(Wisiand = 15ps) a static magnetic island. At this partic¬ 
ular time slices with a static magnetic island. Fig. 5(b), 
5(d), 5(f), 5(h) and 5(j) show the propagation of an ion 
temperature finger-like structure across the island asso¬ 
ciated with a clockwise flow inside the island. This prop¬ 
agation is due to the enhancement of the perpendicular 
heat flux transport close to the X-point and the clockwise 
flow rotation inside the island. At other time slices, we 
can observe opposite characteristics with a radial sym¬ 
metry with respect to x = 0 for lower ion temperature 
structures associated with counter-clockwise flow rota¬ 
tion inside the island. Fig 6 shows the time average of 
the ion temperature over t G [500,3000] a/cg. We ob¬ 
serve the expected flattening of the ion temperature due 
to the effect of the island. However, this flattening is 
not poloidally invariant. The flattening seems to be ex¬ 
tended toward the top of the island. This is due to the 
direction of the poloidal zonal flow which goes along de¬ 
creasing poloidal coordinates at the middle of the box 
(i.e., \x\ < 10) and along increasing poloidal coordinates 
at the edge of the box (i.e., \x\ > 10). Moreover, an ad¬ 
ditional poloidal variation of the flattening is observed in 
Fig. 7 which contains the poloidal average < Tio^Ti >y^t 
of the total ion temperature as well as the profiles across 
the X-point < T^o Ti{y = 6/2) >t and the 0-point 
< Tio ± Ti(p = 0) >t. As a result, the flattening of the 
ion temperature is dominant on the profile across the O- 
point (^ = 0). As explained with the previous figure, this 
is due to the fact that all positive fluctuations of the ion 



xiPs) 


X{ps) 


(a) ^^island — 0 


(b) ^^island — 15ps 





Xips) X{p,) 


(g) Wisiand — 0 


(h) Wisiand — 15ps 



FIG. 5: Snapshots of the total ion temperature Ti = Tio(x) + 
Ti{x,y) at t = {1506,1516,1526,1536,1546}a/cs respectively 
for (a)-(b), (c)-(d), (e)-(f), (g)-(h) and (i)-(j), obtained from 
two simulations with Wisiand = Ops for (a), (c), (e), (g), (i) 
and with Wisiand = 15ps for (b), (d), (f), (h), (j). The same 
colorbar limits (i.e., [0.55; 1.45] and [0.1; 1.9]) are used for all 
time steps of these two simulations (i.e., Wisiand = Ops and 

island — 15ps). 

















































































































































x{ps) ^iPs) 


(a) ^^island — 0 (t>) ^^island — l^Ps 

FIG. 6: Time average on t G [500; 3000]a/cs of the total ion 
temperature Ti = T^o + ft for (a) Wisiand = Ops and (b) 

island — 15ps. 



FIG. 7: Profiles of the total ion temperature Ti = T^o + T. 
The black curve is the poloidal average, the blue (respect, 
green) curve is the prohle across the 0-point at p = Ops (re¬ 
spectively X-point at p = ±32ps) with Wisiand — 15ps. Red 
curve is the poloidal average of the total ion temperature pro¬ 
hle with no island Wisiand = 0. 

temperature (i.e. a hnger-like structure of iou tempera¬ 
ture) which appear close to the X-poiut cuter iuside the 
islaud from the left haud side, follow the maguetic flux 
surfaces toward the right haud side where this Auger dif¬ 
fuses aud cau cross the separatrix toward the right haud 
side. The clockwise rotatiou of the flow iuside the islaud 
is due to the positive radial trausport (toward iucreasiug 
radial positiou) for positive fluctuatious of the iou tem¬ 
perature. lu coutrast, wheu uegative fluctuatious of the 
iou temperature appear close to the X-poiut, the couuter- 
clockwise rotatiou domiuates due to the uegative radial 
trausport (toward decreasiug radial positiou) of uegative 
fluctuatious of the iou temperature. 

Fiually, due to our forciug of the couservatiou ou av¬ 
erage of the backgrouud iou temperature gradieut by 
the preseuce of the viscosity term Vg in the equatious, 
the flatteuiug of the total iou temperature across the O- 
poiut produces au iucrease of the iou temperature gradi¬ 
eut arouud the X-poiut. Moreover, by lookiug at smaller 
variatious of the time averaged total iou temperature pro¬ 


file across the 0-poiut we cau see that arouud the posi- 
tious of the iuuer {dX x = —1bps) aud outer (at x = 15ps) 
separatrix, the gradieuts are iuverted aud become locally 
positive. The effect of these iuversious of the iou temper¬ 
ature gradieut is to geuerate local variatiou of the flow. 
We remark that close to the middle positiou x = Ops, ei¬ 
ther with or without the backgrouud static islaud, a very 
small flatteuiug teudeucy of the iou temperature profile 
appears due to a small maguetic shear. This flatteuiug 
is uegligible iu comparisou to the flatteuiug due to the 
maguetic islaud. 

Iu Fig. 8 we cau see the radial iou heat flux < v^Ti > 
due to the E x B drift averaged over the time rauge 
t G [500; 3000] a/cs for differeut sizes of the islaud 
Wisiand ^ {0;10;20}ps respectively for the figures (a), 
(b) aud (c). With uo static islaud (a), the flux is couceu- 
trated arouud the resouaut surface at x = 0. This time 
averaged radial iou heat flux is almost homogeueous iu 
the poloidal directiou. By iucreasiug the size of the static 
maguetic islaud, we observe a compressiou of the time av¬ 
eraged iou heat flux toward the X-poiut but, more pre¬ 
cisely toward the top of the islaud. ludeed, the poloidal 
profile is more peaked uear y ^ 20ps aud the maximum 
value of this poloidal profile is approximately 3 times 
higher for Wisiand = 20ps thau for Wisiand = Wps aud 
approximately 10 times higher thau for the case with 
uo islaud. We observe that the shape of the radial iou 
heat flux is cousisteut with the previous observatiou of 
the flatteuiug of the iou temperature across the 0-poiut. 
We also remark that the radial iou heat flux average is 
positive because the radial trausport is such that Vx is 
uegative (respectively positive) for uegative (respectively 
positive) iou temperature fluctuatious Ti. To obtaiu more 
detailed uuderstaudiug of the global effects of the islaud 
ou the turbuleuce aud trausport, we ueed to compare the 
radial iou heat flux by averagiug Fig 8, over the poloidal 
dimeusiou y. 

We use a parameterizatiou of the poloidal average (* * * )^ 
of the time-averaged iou heat fluxes of Fig. 8. The pa¬ 
rameterizatiou is doue usiug a flttiug of the profiles with 
the followiug Gaussiau 

fix) = fo + fi exp 

The variatiou of the parameter /o describes the effects of 
the islaud outside the islaud (i.e., \x\ > Wisiand)^ fo + fi 
the effects arouud the positiou x = 0 iuside the islaud, 
aud Wf the characteristic width of the Gaussiau pro¬ 
file. Fig. 9 represeuts the variatious of these parameters 
used by the Gaussiau flttiug of the radial iou heat flux 
Q = VxTi correspoudiug to the trausport property aud 
the squared electric poteutial correspoudiug to the 
turbuleuce property. 

We observe the fact that the width of the Gaussiau pro¬ 
file fits as well as the maximum value at x = 0 iucrease 
for the trausport aud the turbuleuce beyoud a threshold 
islaud size. The threshold is arouud bps for the traus¬ 
port aud betweeu bps aud lOps for the turbuleuce. These 
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^(Ps) 


(a) Radial heat flux for Wisiand = OPs 



^(Ps) 

(b) Radial heat flux for Wisiand = lOps 


xlO“^ 



-20 0 20 


(c) Radial heat flux for Wisiand = 20ps 

FIG. 8: Comparison of the time average of the radial heat 
flux Q = VxTi over the time range t G [500,3000] a/cs 
between the following sizes of the static magnetic island 
Wisiand G {0;10;20}ps. The white plain (resp. dashed) 
curves are the island separatrix (resp. 15 iso-contour mag¬ 
netic flux surfaces). The radial ion heat flux is peaked inside 
the island between the X-point and the 0-point. 




(a) Width Wq of the fitting 
Gaussian profile 


(b) Max Qo QI and min Qq 
of the fitting Gaussian profile 



°0 5 10 15 20 

^^islanJyPs ) 

(c) Width W(^ of the fitting 
Gaussian profile 


0.08 



^^islandLPs ) 


(d) Max 00 + 01 and min 0o of 
the fitting Gaussian profile 


FIG. 9: Comparison of the Gaussian fitting parameters 
against the size of the static magnetic island Wisiand G 
{0; 1; 2; 5; 10; 15; 20}ps for (a)-(c) the radial ion heat flux and 
for (b)-(d) the square of the electrostatic potential. The time 
range is t G [500,3000] a/cs. Beyond a threshold around 
Wisiand ~ 5ps, the static island globally enhances the trans¬ 
port and turbulence even outside the island. 


thresholds are about the same order as the characteristic 
length of the turbulence bps to lOps of our simulations. 
Moreover, the parameters Qo and correspond to the 
effect of the island respectively on the transport and the 
turbulence outside the island. Even if it is a small vari¬ 
ation, we observe an increase of the transport and the 
turbulence far away from the island. This observation 
means that there is an enhancement of the turbulence 
and the transport independent of a simple displacement 
by the presence of the island. 

Finally, these results given by a simple poloidal average 
can be verified by using another average which conserves 
the structure of the static magnetic island. To meet this 
goal, a last investigation is to try to recover the pre¬ 
vious results by using the flux surface average (• * * )^^ 
instead of the poloidal average (• • • )^. Fig. 10 represents 
the flux surface averages of the time-averaged radial ion 
heat flux for different sizes of the magnetic island. The 
first figure corresponds to profiles of the radial ion heat 
flux inside the island versus the flux surface coordinate 
— pdo)/{'ipsep — for different sizes of the magnetic 
island and the second figure to profiles outside the island 
versus {ip - 'lpsep)/{'lpa - i^sep)‘ The values {'Ipo.'lpsep.'^a) 
of the magnetic flux correspond respectively to the O- 
point, the separatrix, and the last closed flux surface in 
our simulation box. Beyond the last closed flux surface, 
all integrals along each opened flux surface linearly de¬ 
crease to 0, which is their minimum values at the corner 
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of the simulation box of the flux surface corresponding 
to '0eg(—As a summary, for Fig. 10(a) (respec¬ 
tively Fig. 10(b)) corresponding to the radial ion heat 
flux inside (respectively outside) the island, the abscissa 
coordinate represents the 0-point (respectively the sepa- 
ratrix) at the value 0 and the separatrix (respectively the 
last closed flux surface) at the value 1. We observe inside 




(b) Flux surface average of the heat flux outside the island. 


FIG. 10: Comparison of the time and flux surface averages 
of the radial ion heat flux for the following sizes of the static 
magnetic island Wisiand ^ {5; 10; 15; 20}ps. The effect of the 
poloidal asymmetry inside the island and the threshold of 
Wisiand ~ 5ps are consistent with these flux surface averages. 

the island (Fig. 10(a)) an increase of the shift of the ra¬ 
dial ion heat flux toward the separatrix instead of being 
homogeneous everywhere inside the island. However, the 
asymmetry between the top and the bottom of the island 
cannot be recovered in comparison to what have been 
shown above with the 2D slices. With the observation of 
the radial ion heat flux outside the island (Fig. 10(b)), the 
background transport increases beyond a threshold size 
of the magnetic island between bps and lOps, related to 


the characteristic length of the transport with no island. 
This enhancement of the transport outside the island is 
due to the filaments which appear at different time slices 
as described previously. As an example with a static 
magnetic island Wisiand = 20ps, we typically observe fil¬ 
aments of the positive fluctuations of the ion temperature 
starting around the position (x,p) = (—10,25), crossing 
the separatrix, following the magnetic flux surfaces until 
the position (x, y) = (20, 0) and finally diffusing outside 
the island in the region of lower ion temperature. 



x{Ps) x{pP 

(a) {b>^ for Wisiand = Ops. (b) {b>^ for Wisiand = Wps- 


FIG. 11: Gomparison of the time averaged electric potential 
(0)^ on t G [500; 3000] a/cs for a static island of width (a) 

W^island — Ops and (b) Wisiand — 15ps. 



(a) Mode ky = 0 for 
W^island — OPs- 


(b) Mode ky = ^ for 
W^island — 15ps. 



(c) Mode ky = ±27r/6 for 
W^island — OPs- 


(d) Mode ky = ±27r/6 for 
W^island — 15ps- 


FIG. 12: Reconstruction of the 2D electric potential from one 
selected poloidal mode: (a)-(b) from ky — b and (c)-(d) from 
ky = ±27r/5 with (a)-(c) no island and (b)-(d) Wisiand = 15ps. 


It is important to note that in addition to the turbu¬ 
lent heat flux, there is also transport due to static E x B 
flows which arise due to the presence of the island, that 
will also drive transport. Fig. 11 shows the time-averaged 
electrostatic potential for the cases of Wisiand = Ops and 
15ps, respectively. In the no-island case, the mean po¬ 
tential corresponds to a ky = 0 zonal flow with radial 
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wavenumber approximately equal to nja. However, in 
the finite island case (Fig. 11(b)), the mean potential 
has both significant ky = {) and ky = 27r/b components, 
with comparable magnitude to the no-island case. These 
components are separately visualized in the (x,y) plane in 
Fig. 12. Note in particular that the finite ky equilibrium 
flow in the finite island case is not symmetric about the 
0-point at ^ = 0, and thus contributes to the asymmetry 
in radial heat flux shown in Fig. 8. 

As a conclusion of our investigation of the effects of a 
static island on turbulence and transport, we observed 
the flattening of the ion temperature profile across the 
0-point {y = 0) due to the island as well as an increase 
of its gradient across the X-point. The increase of the 
ion temperature gradient around the X-point enhances 
the ITG instability and drives more turbulence in this 
region. Due to the radial transport, positive or negative 
structures of the ion temperature enter inside the island 
and flow along the magnetic flux surfaces. Inside the is¬ 
land, the time average of the radial ion heat flux is then 
peaked on the top of the island due to the fact that the 
flow naturally goes along decreasing poloidal coordinates 
around the radial position x = 0 (as observed with no 
island). This privileged direction of the flow is due to 
our choice of the sign of the equilibrium poloidal mag¬ 
netic flux (i.e., the sign of the magnetic shear). Finally, 
because of the presence of the island, the turbulence and 
transport are significantly increased when the size of the 
static magnetic island is bigger than the characteristic 
size of the turbulence and transport. 

To make progress toward understanding the self- 
consistent interaction, some results of the feedback of 
the modified turbulence on the dynamical evolution of 
the island is investigated below. 


IV. FIRST RESULTS TOWARD THE 
FEEDBACK ON THE ISLAND 

One of our interests is to understand the self-consistent 
interaction between fast dynamics of microturbulence 
and the slow dynamics of a magnetic island. In the pre¬ 
vious section we focused on the effect of a static magnetic 
island on transport and ITG microturbulence. From 
these results, we report here a quantitative study of the 
effect of microturbulence on the slow dynamics of the 
magnetic island. 


A. Formulation of a turbulent resistivity 


We begin by adopting a mean-field type approach to 
exploit the separation of time and spatial scales between 
the turbulence and the island. We can then rewrite the 
magnetic flux evolution equation as 


Pdt^i 


By^ody(j)i ydjzj T ^ 


— n 


, (30) 


where By^i^dy = Vy = Vy — Vy, the second term is the 
usual classical resistive current and the last term corre¬ 
sponds to the average of the Ohm’s law stress. This last 
term can be written as a coherent part —rjtnvhjzj and 
a non-coherent part —t^ncJnc with respect to the island 
current where the effective turbulent resistivity describes 
the self-consistent feedback of the turbulence on the dy¬ 
namics of the island 


'^turb 


JzJ 


x,y 





(31) 


Moreover, in addition to the dynamical evolution of the 
magnetic flux of the island a dynamical evolution 
of the electric potential associated with the magnetic 
island has also been developed. Here, we do not focus on 
this part because in our previous simulations we imposed 
only a poloidal magnetic flux and we omit the presence 
of an associated electric potential. Indeed, the associated 
electric potential of the island is a very slow dynamical 
quantity in comparison to the fast evolution of (j) and 0/ 
has almost no effect on the dynamical evolution of the 
magnetic island poloidal flux. 


B. Evaluation of turbulent resistivity in simulation 


For the dynamical equation of the island, simulations 
can be challenging since we need to resolve the fast 
time scales of the microturbulence as well as the slow 
time scales of the evolution of the magnetic island and 
the tearing mode. In this work, we have exploited the 
timescale separation between the island and turbulence 
dynamics to approximate the island as fixed, consistent 
with previous studies by other groups. However, we can 
still begin to probe how we might expect the turbulence 
to act on the island in a fully dynamic model by calculat¬ 
ing the time-averaged forcings of the turbulence on the 
imposed island structure. 


Fig. 13 represents the Ohm’s Law stresses 


t/;, 0 — h 


averaged over t G [500; 3000] a/cg for different sizes of the 
static magnetic island Wisiand ^ {1; 2; 5; 10; 15; 20}ps. 
There is an obvious link between the structures of these 
Ohm’s Law stresses and those of the radial ion heat fluxes 
shown in Fig. 8. The effective turbulent resistivity r/turb 
given by the Eq. (31) is then computed from the extrac¬ 
tion of the coherent part of the Ohm’s law stress with re¬ 
spect to the static current associated with the magnetic 
island. The ratio between the effective turbulent resis¬ 
tivity and the classic resistivity rjcX is shown in Fig. 14. 
The first result is the expected fact that the effective 
turbulent resistivity, which quantifies the effects of the 
feedback of the turbulence on the dynamical evolution of 
the magnetic island, decreases by 1 order of magnitude 
as the island size increases. The second result is the am¬ 
plitude scale between the effective turbulent resistivity 
and the classic resistivity. Indeed, for all island sizes, the 
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FIG. 13: Ohm’s law stresses 0 — hj averaged over t G 

[2250; 3000] a/cg for different sizes of the static magnetic is¬ 
land Wisiand G {1; 2; 5; 10; 15; 20}ps (from left to right and top 
to bottom). The black plain (resp. dashed) curves are the is¬ 
land separatrix (resp. 15 iso-contour magnetic flux surfaces). 
The structures of the Ohm’s law stresses are consistent with 
the radial ion heat fluxes. 



FIG. 14: Ratio of the effective turbulent resistivity ?7turb given 
by Eq.(31) over the classic resistivity rjci versus the size of the 
island Wisiand- 


effective turbulent resistivity is at least 1 order of mag¬ 
nitude higher than the classic resistivity. However, even 
with a difference of 2 orders of magnitude for small is¬ 
lands, the total effect of the self-consistent feedback of 
the turbulence on the island is given by the multiplica¬ 
tion of the effective turbulent resistivity 77turb with the 


current generated by the magnetic island jzj- The later 
is approximately 10 times smaller for small island sizes 
than for larger island sizes. This investigation clarifies 
some expected mechanisms toward the self-consistent in¬ 
teraction between the ITG microturbulence and the dy¬ 
namical evolution of a magnetic island. Perspectives are 
discussed below. 


V. DISCUSSIONS AND FUTURE WORK 

In this work, we have studied the dynamics and asso¬ 
ciated thermal transport due to ITG microturbulence in 
the presence of static magnetic islands of varying width. 
We find that the magnetic island is responsible for en¬ 
hancement of the turbulent fluctuation amplitudes and 
radial heat flux, which increases rapidly above a thresh¬ 
old island width comparable to the correlation length of 
the turbulence in the no-island limit. Broadly consistent 
with theoretical expectations, we observe a flattening of 
ion temperature inside the island, and localization of the 
radial heat flux near (but not symmetric about) the is¬ 
land X-point. It should be noted that this flattening 
occurs due to turbulent transport processes and in-plane 
flow along flux surfaces, as the simulations presented here 
do include explicit conductive parallel heat fluxes of the 
form II = —/c||V||T^ that are often used to represent 
both parallel collisional dynamics and Landau-damping 
effects. Thus, the combination of self-consistently calcu¬ 
lated radial and poloidal turbulent heat fluxes in the pres¬ 
ence of the imposed island (along with the self-consistent 
electrostatic flow which also forms) are sufficient to lead 
to temperature flattening and equilibration on the total 
magnetic flux surfaces. A novel observation has been 
made that a asymmetric shift in location of the peak ra¬ 
dial heat flux relative to the X-point is due to inclusion 
of self-consistent turbulence effects (such as the inherent 
phase velocity of the turbulence) which are not included 
in theoretical studies which simply approximate the tur¬ 
bulent fluxes in terms of anomalous diffusion coefficients. 
In addition to calculating the response of the turbulence 
to the presence of the island, we also investigated the 
back-reaction of the turbulence on the island, using a 
mean-field approach to calculate an effective turbulent 
resistivity which would act on the island structure. For 
the parameters considered, our simulations predict that 
while this turbulent resistivity is always at least an or¬ 
der of magnitude larger than the collisional value used, 
it is most effective at small island widths and decreases 
quickly with increasing island size. 

There are a wide variety of future directions to be pur¬ 
sued in future work, building upon these results. Fore¬ 
most are inclusion of parallel heat flux due to both col¬ 
lisions and Landau damping, as well as fully 3D effects, 
to make closer connection with experiments. Equally im¬ 
portant in our opinion is to transition from statically im¬ 
posed islands to self-consistently evolving islands driven 
by tearing unstable current profiles. An important nu- 
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merical challenge to be solved for these simulations in 
3D is the implementation of preconditoners and implicit 
advance schemes suitable for long-time integration of the 
slowly evolving island in the presence of fast ITG turbu¬ 
lence and even faster damped Alfven waves in 3D geom¬ 
etry. Greater understanding the self-consistent transport 
enhancements and profile flattening effects, as well as tur¬ 
bulent forcings of the island are needed to develop a fully 
predictive model of coupled tearing mode and turbulence 
dynamics. We can then begin to include more realistic 
geometric effects in the equilibrium structure, to prop¬ 
erly describe the coupled turbulence and reconnection 


processes which drive observed magnetic islands. 
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